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Abstract 



We study the pecuhar motion of non-relativistic matter in a fuUy covariant way. The exact 
nonhnear equations are derived and then apphed to the case of pressure-free matter, moving 
' relatively to a quasi-Newtonian Eulerian frame. Our two-frame formalism facilitates the 

study of the nonlinear kinematics of the matter, as the latter decouples from the background 
expansion and starts to "turn around" and collapse. Applied to second perturbative order, 
our equations provide a fully covariant formulation of the Zeldovich approximation, which 
^\ • by construction addresses the mildly nonlinear regime of structure formation. Employing a 

, dynamical system approach, we show that, just like in the Newtonian case, the relativistic 

' treatment also predicts that pancakes are the natural end-structures for any generic over- 

, density. 

P-t', PACS number(s): 04.25.Nx, 98.80.Hw 

O 

c/3 '. 1 Introduction 

Studies of linear perturbations are crucial for understanding the way structure formation has 
^ , progressed in our universe. In the standard model density fluctuations grow slowly, via gravita- 

' tional instability, to form the galaxies, the clusters of galaxies, the superclusters, the filaments 

^ and the voids that we see in the universe today. Despite the simplicity of this idea, however, our 

analytical understanding of the advanced stages of gravitational collapse is still limited. The 
reason is that the linear theory is a good approximation only at the initial stages of the collapse, 
when the density contrast of the perturbation is well below unity. Most of the structures in the 
universe, however, have density contrasts well in excess of unity. To understand the evolution of 
these objects we need to go beyond the limits of the linear approximation. The inherent com- 
plications of nonlinear collapse, however, mean that analytical studies are possible only when 
certain simplifying assumptions are made. An approximate approach that extends well into the 
nonlinear regime, up to the virialization of bound objects, is the spherical collapse model [1] 
(see also [2]). Although the spherical model became popular because of its simplicity, in reality 
it stops short of explaining key features of the observed universe. Gravitational collapse does 
not seem to proceed isotropically. All galaxy surveys show structures with complicated triaxial 
shapes, which require a nonspherical analysis if they were to be explained. 
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The Zeldovich approximation is not restricted to spherical symmetry [3] . It is a kinematical, 
Lagrangian approach that addresses the issue of anisotropic collapse by extrapolating into the 
nonlinear regime a well known linear result (see [4] for a review and also [5] for related discus- 
sion) . More specifically, the Zeldovich ansatz assumes that the acceleration- free and irrotational 
linear peculiar motion of pressure-free matter also holds during the nonlinear collapse. The con- 
sequences of this hypothesis are quite dramatic. What Zeldovich showed was that any generic 
(i.e. non-spherical) ovcrdensity will undergo a phase of anisotropic, effectively one-dimensional, 
collapse leading to the formation of two-dimensional flattened structures that are widely known 
as "pancakes". Over the years, the Zeldovich approximation has provided a great deal of in- 
sight into the initial nonlinear evolution of density fluctuations. It has also inspired a number 
of developments concerning the behavior of the velocity field in irrotational flows (see [6] for 
a representative list).^ Most of the available studies, however, are Newtonian and, although 
there are general relativistic approaches in the literature (see [12]), a fully covariant treatment 
is still missing. Moreover, Newtonian and relativistic collapse seem to move further apart in the 
presence of anisotropy (e.g. see [13]). In the Newtonian approach pancakes are the natural at- 
tractors, at least in the Zeldovich approximation; a result that seems well confirmed by N-body 
simulations [14]. The relativistic "silent universe" models, on the other hand, point towards 
spindle-like singularities [15, 16].^ In other words, the final fate of a collapsing overdensity re- 
mains an open issue. In the present article we attempt to address this question by looking at 
the dynamics of nonlinear collapse in a perturbed Einstein-de Sitter universe. 

One issue that emerges naturally, especially when attempting a relativistic analysis of pe- 
culiar motions, is that the associated velocities should be defined relative to a preferred frame 
of reference. This frame is necessarily non-comoving with the matter, since by definition there 
are no peculiar velocities relative to the matter frame, which is the truly Lagrangian frame. It 
should be emphasized that in the standard literature the term "comoving" refers usually (and 
sometimes misleadingly) to a fictitious background 4-velocity, rather than to the actual velocity 
of the matter. In the covariant approach one seeks to define the aforementioned preferred frame 
in a physical way.^ Motivated by [20]- [22], we chose a 4-velocity field that is both irrotational 
and shear-free, thus defining a quasi-Newtonian frame that closely corresponds to Bardeen's 
quasi-Newtonian gauge [23]. Our aim is to set up the general framework for the fully covariant 
and fully nonlinear treatment of non-relativistic peculiar motions. We then exploit the natural 
transparency of the covariant formalism to provide a fully covariant version of the Zeldovich 
approximation. As is well known, the latter addresses the early stages of nonlinear collapse. We 
call this period the "mildly nonlinear regime" . Given that we deal with the post-recombination 
era, we assume an Einstein-de Sitter background and consider the second order evolution of 
the perturbed variables. In the frame of the pressure-free matter, we find that the equations 
describing the collapse of the non-relativistic component to second perturbative order do not 
include any tidal effects. They reduce to a set of ordinary second order differential equations, 

^The limits of the Zeldovich approximation, both linear and nonlinear, have been investigated in [7]. Some 
authors have suggested modifications of the standard Lagrangian approach to increase its accuracy, by adding 
vorticity [8], or viscosity [9], while others have treated either the velocity potential as a constant [10], or the 
gravitational potential itself [11]. 

^Silent universes are inhomogeneous spacetimes, with irrotational pressure- free matter and zero magnetic Weyl 
tensor. By construction, they cannot support any communication through sound or gravitational wave signals, 
which explains the term silent. This class of models contains the well known Szekercs solution as a special case. 
For an introduction and further discussion on silent universes we refer the reader to [17]. 

^An extensive presentation of the covariant approach to cosmology can be found in [18]. For an updated review 
the reader is referred to [19]. 
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which can be rescaled into a planar dynamical system with one-dimensional pancakes as its 
natural attr actors. 

2 The non-linear peculiar kinematics 

2.1 The peculiar velocity field 

Velocity perturbations and their implications for gravitational collapse are fundamental in any 
structure formation scenario. In the post-recombination era and on scales well below the Hubble 
length, the Newtonian theory can adequately address the question of gravitational instability. 
As we are successively probing significant fractions of the Hubble radius, however, the need for 
a relativistic treatment grows. In what follows we will lay down the formalism for the fully 
covariant treatment of peculiar velocities in cosmology. 

To begin with, consider a family of timelike worldlines tangent to the reference 4-velocity 
vector field Ua, normalized so that Uau"' = — 1, and assume the presence of matter moving with 
4-velocity Ua so that 

Ua = l{Ua + Va) ■ (1) 

Here Va is the peculiar velocity of the matter with respect to Ua, 7 = (1 — is the Lorentz- 

boost factor between the two frames and = VaV"'. Note that Va vanishes in the background, 
because in a Robertson- Walker model both Ua and Ua would be chosen parallel to the canonical 
time direction, which makes it a first order gauge-invariant quantity. Also, for non-relativistic 
peculiar motions, like those we will be dealing with, v'^ <^ 1. This means that 7 2± 1 , which in 
turn ensures that UaU°- 2± — 1. Apart from this restriction, the frames Ua and Ua are completely 
general. Later, however, they will be respectively identified with the quasi-Newtonian observers 
and with the (pressure- free) matter (see Sec. 3). 

The 4-velocity field Ua of the matter can be used to define an additional time direction and 
an associated projection tensor given by 

Kb = 9ab + UaUb (2) 

where Qat is the spacetime metric. Note that Va is not orthogonal to Ua- Indeed, even for 7=1, 
Eqs. (1) and (2) ensure that Uav"' = ^ Q and that habV^ = Va + v^Ua 7^ Va respectively. 

2.2 Decomposition of the peculiar velocity gradients 

Following the standard covariant splitting of the 4-velocity gradients into the irreducible kine- 
matical quantities of the motion, we will also decompose the gradients of the peculiar velocity 
field. We take Ua as our time direction because this leads to simpler equations later on (in effect 
we are choosing a non hypersurface orthogonal threading of the space-time manifold). We then 
arrive at the expression 

VbVa = BbVa - VaUb + Uai^bc - ^bc + ^Qhc)'"" - UgS b^^ - {v^)'UaUb , (3) 

where Va = u^^^^a and = hJ'Vb is the covariant derivative operator orthogonal to Ua- It 
should be emphasized that 0, dab ^ab a^re the expansion scalar, the shear and the vorticity 
associated with Ua- Also, in deriving Eq. (3) we have used the relation Uav'^ = v^, which holds 
for 7 = 1. 
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Similarly, the projected gradient T)i,Va of the peculiar velocity is decomposed into the irre- 
ducible quantities of the motion as follows 

f)bVa = ^Qhab + ^ab + <^ab , (4) 

where^ 

e = Da?;" , Gab = t>{bVa) and Uab = t>[bVa] ■ (5) 

The above respectively represent what one might call the expansion (or contraction) scalar, the 
shear tensor and the vorticity tensor associated with the peculiar motion. Clearly, when 6 > 
we are dealing with an expanding fluid element, whereas B < indicates collapse. 

2.3 The peculiar motion 

The motion of the matter is monitored through a set of two fully nonlinear propagation equations, 
with time derivatives taken along the matter flow lines. The first determines the evolution of 
the peculiar velocity vector 

V{a) = Aa- A(a> - - ^)v{a) " {^ab " Crab)v^ - {oJab " u;ab)v'' , (6) 

and the second governs the behavior of its projected gradients 

ha%'^ (OrfWc) = BbVa - + OJ^b + ^Qh^b^ (aca + iOca + ^©^ca) + ^6 

-Aa (acb + OJcb + \Qhcb) + AJ^bV^ - K'^h'^RecsdV^u' . (7) 

The former of the two equations can also be seen as the transformation law of the 4-acceleration 
between the "tilded" and the "non-tilded" frame, with Aa = vf'Vb^a and Aa = u^Vb^a- The lat- 
ter comes after applying the Ricci identity to Va-, projecting orthogonal to Ua and then employing 
decompositions (3) and (4). 

The trace, the symmetric trace-free part and the antisymmetric component of Eq. (7) provide 
the nonlinear evolution equations for the irreducible kinematical variables of the peculiar motion. 
In particular we have 

- [^ab - ^ab + l&Kb) A^' + i"D„i;2 (8) 

for the peculiar expansion (or contraction), 

h{a'^hb/^cd = -5(0^06 + ©CTab) - ac{ahf - ojc{a'^b}'^ + D(a?)fe) - A(^ahb}'^Vc - h^a'^hb/RgcsdV^U^ 
-^cia'^b)" - ^da^b)" - Al^a {^b)c " '^6>c + 50^6>c) + ^(aDfo)?;^ (9) 

for the peculiar shear and 

hla^hbfu^cd = -l{&<^ab + Qi^ab) + ^c[a^bf + ^c[a^bf " '^[ah] " A^ahbfVc - ^a^hf Recsdv"" 

+^c[a^b]' + ^ciahf - \ (o-fejc " ^b]c + g0^6]c) v'' + i[aDb]^^^ • (10) 

for the peculiar vorticity. 

*Round brackets indicate symmetrization and square ones antisymmetrization. Angled brackets, on the other 
hand, indicate the projected, symmetric and trace- free part of tensors and also the projected component of vectors. 
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2.4 The spacetime curvature effects 

The curvature terms in the right hand side of Eqs. (8)-(10) convey the gravitational effects on 
the peculiar kinematics. To evaluate their contribution we consider the standard decomposition 
of the Riemann tensor 

Rabcd = Cabcd + \{gacRbd + QbdRac — 9bcRad — 9adRbc) — ^R{9ac9bd — 9ad9bc) , (H) 

where Rab = R^acb R = R^^a are the Ricci tensor and the Ricci scalar respectively, while 
Cabcd is the Weyl tensor. The latter describes the long range gravitational field, namely tidal 
forces and gravity waves. Relative to Ua the Weyl tensor decomposes further as follows [24] 

Cabcd = UaUcEhd — UaUdEijc — UbUcEad + UbUdEac — Eabe^cdfE^^ 

+eabeUcH''d - SabeUdH^'c + EcdfUaHb''' - EcdfUbHa^ , (12) 

where Eab and Hab are the associated electric and magnetic Weyl components. Note that 
^abc = Vabcd'u'^ is the totally antisymmetric 3-dimensional permutation tensor orthogonal to Ua, 
with r]abcd being the spacetime alternating tensor. Introducing the above into expression (11) 
we obtain 

RabV^vl" = -QaV" - [^(/X + 3p) - h\ v\ (13) 
h{a%)'^RecsdV^u' = -v'^Eab + £cd(aV^H'^b) + ^T^ab + Uiahf^c , (14) 

and 

h[a%fRecsdV^u' = e^d[aV''H'^b] + k^laKfVc , (15) 

where /i, p, Qa and tTqj are respectively the energy density, the isotropic pressure, the heat flux 
and the anisotropic stresses of any matter fields that may be present as measured in the Ua 
frame. 

2.5 The fully nonlinecir equations 

Substituting results (13)-(15) into Eqs. (8)-(10) we obtain the following fully nonlinear expres- 
sions for the irreducible kinematical variables 

^ = -^ee~daba''^+C[^ab^'''' + BaV'' + A''Va + qaV''+[^iil + Sp)-A]v^ 

- {^ab - ^ab + l&hab) A'^v'' + A^DaV^ , (16) 

h(a''hb/^cd = -ii&^ab + Q^ab) - ^cia^f " ^c{a'^b)'' + '^{a'"b) " Ai^ahf-^c + V^Egb - S cdia^" H"^ b) 
-^V^T^ab - \,q{ahb)''Vc - ^c{a'^b)'' " ^ciah)" " ^(o {^b)c " '^b)c + ^Q^c) """^ 

+Ai^a^b)v''. (17) 

and 

hla'^hbf'^cd = -\{Q<^ab + Qi^ab) + 0-c[a^bf + i^c[a<^bf " i^[aVb] " AiahfVc - S cdia^" H"^ b\ " iQiohbfVc 
+^c[a^b] + i^c[ahf - \a (c^6]c " ^b\c + l^h]^ + AyJ)b\V^ . (18) 
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The above, together with the propagation equation of Va itself, that is 



V(a) = Aa- - 5(6) " 0)w(a> " i^ab " '^ab)^'^ " {l^ab " l^ab^ 



,b 



(19) 



determine the pecuUar motion completely. Finally, we remind the reader that to this point we 
have not explicitly identified Ua and Ua with any physical frame. In other words, the generality 
of formulae (16)-(19) is only restricted by the fact that the relative motion of the aforementioned 
observers is non-relativistic (i.e. 7 = 1). 

3 Application to pressure-free matter 

3.1 The quasi-Newtonian observers 

When it comes to studies of relativistic cosmological perturbations, Bardeen's gauge-invariant 
theory is the most popular approach [23]. Various gauge choices are possible in this formalism. 
The longitudinal gauge for example, which assumes constant-time slices with zero-shear normals, 
sets up a frame of reference that emulates the Eulerian observers of the Newtonian theory, thus 
motivating the term "quasi-Newtonian" gauge. Since this frame is non-comoving, one can use 
it to study the relative velocity effects. The covariant alternative mimics the fixing conditions 
of Bardeen's quasi-Newtonian gauge by introducing an l-|-3 threading of the spacctimc via 
a 4- velocity field that is both irrotational and shear- free [21, 22]. Thus, covariantly, quasi- 
Newtonian cosmologies are almost-FRW universes equipped with a family of observers, moving 
non-relativistically relative to the comoving observers, whose motion is also non-rotating and 
shcar-frcc. In addition to, such cosmological models do not support gravitational waves, which 
further justifies the term quasi-Newtonian [21]. 

In the current section, the quasi-Newtonian observers will be identified with the 4-velocity 
vector field Ua defined in Sec. 2. As a result, the shear and the vorticity tensors associated with Ua 
will vanish (i.e. (Tab = = i^ab)- The rest of the kinematic, dynamic and gravito-clcctromagnctic 
quantities depends on the choice of the matter component and on the transformation laws 
between the Ua and the Ua frames. For the full list of the exact nonlinear relations the reader is 
referred to [25, 26]. 

3.2 The nonhnear pecuhar motion of the matter 

So far, Va is an arbitrary peculiar velocity field which has not been associated with any par- 
ticular matter source. To this point, Eqs. (16)-(19) are completely general and 7 = 1 is the 
only restriction imposed. Next, we will align with the motion of a non-relativistic mat- 
ter component. This choice allow us to address the question of velocity perturbations in the 
post-recombination era, when the universe is dominated by pressure-free dust, while it con- 
siderably simplifies Eqs. (16)-(19). In particular, = in the matter frame by definition 
and p = TTab = ^0 = for dust. Therefore, the nonlinear peculiar motion of the pressureless 
component is governed by 



V{a) = -A{a) - - ^)v{a) - {^ab - ^ab)v^ - {i^ab - <^ab)v^ , 



(20) 



e 



(21) 



^ab = -\{^^ab + ^CFab) - ^da^b)" " i^cia^b)" + ^{aVb) + Kb - E cd{aV'' H'^ b) 



(22) 
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and 



+^c[a<^b] +'^c[ah] ■ (23) 

Note the effect of the long range gravitational field on the anisotropy of the collapse, which 
propagates through the electric and magnetic Weyl terms in Eq. (22). Their coupling to the 
peculiar velocity of the matter means that the tidal effects are of higher pcrturbative order than 
the rest. This will prove crucial when Eqs. (20)-(23) are applied to the mildly nonlinear regime 
of gravitational collapse (see Sec. 4.1). Also note how the tidal field acts as a source of peculiar 
vorticity through the magnetic Weyl component in Eq. (23). 

3.3 The linear peculiar motion of the matter 

For quasi-Newtonian observers Gab = = Uab and the linear transformation laws between the 
quasi- Newtonian and the matter frames guarantee that [25, 26] 

^ab = ^ab , i^ab = ^^ab and 6 = 6 + 6. (24) 

Also to first order, the matter density and the gravito-electromagnetic quantities are related by 

/i = , Eah = Eab and = . (25) 

Finally, the 4-acceleration Aa of the quasi-Newtonian observers takes the form [22] 

^a = Da0, (26) 

where (j) is the peculiar gravitational potential. The latter has the trivial gauge freedom to add, 
or subtract, an arbitrary homogeneous function. Note that D^^ = D^^ to first order, given 
that hab = hab to zero order. Employing the above results Eq. (20) reduces to the following 
expression 

Va + iQVa = -Dafl) , (27) 

which monitors the linear evolution of the peculiar velocity between the quasi-Newtonian and 
the matter frames relative to the latter. Recall that an overdot indicates covariant differentiation 
along Ua (i.e. Va = U^VtVa)- 

For an Einstein-de Sitter background and relative to the quasi- Newtonian observers, the 
linear peculiar velocity has a single growing mode given by v"' = v^a^^'^ [22].^ Here, is the 
peculiar velocity at some initial time and the scale factor has been normalized so that gq = 1. 
Now, to linear order, it is straightforward to employ Eq. (1) and show that u'^'VbVa = u^VbVa, 
which immediately implies that Va has the same linear evolution (i.e. Va oc a^/^) in both the 
quasi-Newtonian and the matter frames. On these grounds, it helps to define the rescaled 
peculiar velocity field 

Va = a-^l\a. (28) 

®The linear evolution law Va oc a}^'^ obtained in [22] assumes that the peculiar velocity is aligned with the 
peculiar acceleration, namely that Va oc Aa. This assumption, which is also part of the Zeldovich approximation, 
ignores contributions from any component of Va that happens to be orthogonal to Aa. As has been pointed out 

in [22], this is not unreasonable since: (i) if the orthogonal component is zero at some initial time, it remains 
so; (ii) if there is an initial orthogonal component, then the expansion and the matter aggregation will serve to 
decrease it. 
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The advantage of is that, to leading order, it is constant both in the quasi-Newtonian and 

in the matter frame by construction. Thus, Va = to hnear order. An immediate consequence 
is that the rescalcd peculiar velocity derives from a potential. Indeed, by substituting (28) into 
Eq. (27) we arrive at the linear expression 

Fa + ieFa = -ieD„$, (29) 

where $ = 20/a^/^0 is the rescaled peculiar gravitational potential. The above clearly manifests 
that 

Va = ^ F„ = -D„$, (30) 

namely that to first order the rescaled peculiar velocity is simply the gradient of the rescaled 
peculiar gravitational potential. Furthermore, in the absence of background rotation, result (30) 
also guarantees that 

D[6K] = -D[feD,]$ = 0, (31) 

to leading order. Therefore, one can always treat the linear peculiar motion of the dust compo- 
nent as acceleration-free and irrotational. 



4 The second order peculiar motion of the matter 
4.1 The original kinematical equations 

The linear equations of the previous section are adequate in the early stages of gravitational 
aggregation when the perturbed quantities are still relatively small. As the perturbations grow 
stronger, however, the linear approximation brakes down and one has to incorporate nonlinear 
effects. During the transition from the linear to the fully nonlinear collapse, a period which one 
might call "the mildly nonlinear era", one can adequately monitor the perturbed variables by 
employing the second-order equations instead of the fully nonlinear ones. In what follows, we 
will employ these second order equations to study the growth of velocity perturbations during 
the early nonlinear stages, as the inhomogeneity decouples from the background expansion, and 
starts to "turn around" and collapse under its own gravity. 

Maintaining the Einstein-de Sitter background of the previous section, we find that the 
second-order peculiar motion of the pressure-free matter is monitored by the following set of 
equations^ 

V(a} = -A(a} - |0^^(a> , (32) 

e = -ie2-2a2 + 2w2^D„{;" + i/x?;2, (33) 

(^ab = -l&'^ab + '^d-c[aUJbf -B[aVb] ■ (35) 

These are obtained from Eqs. (20)- (23) on using the linear transformation laws (24) between the 
kinematical quantities of the quasi-Newtonian and the matter frame. We should also emphasize 

®We employ the linearization scheme used in [26]. In brief, terms that vanish in the background, such as the 
peculiar velocity or the electric Weyl tensor for example, have perturbative order 0{c). where e is the smallness 
parameter. To incorporate second order effects, we include terms of order C(e^, ev, v^) and neglect those of order 
0(e^, ev^, v^) and higher. For more details the reader is referred to [26]. 
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that the system (32)-(35) applies to subhorizon scales only, where the peculiar motion dominates 
over the background expansion (i.e. for O <C 0). 

Note that, to second perturbative order, the peculiar shear evolves unaffected by tidal forces, 
since the electric and magnetic Weyl tensors no longer contribute to Eq. (34). The absence of 
a shear-tide coupling is a direct consequence of the two-frame approach we have employed. It 
means that the evolution of the fluid element proceeds unaffected by the long range gravitational 
field. This in turn has crucial implications for the asymptotic final shape of the collapsing 
overdensity. As we shall see later, the absence of the electric tidal field means that generic 
collapse does not evolve towards the spindle-like singularity that is usually associated with the 
silent universes models [15, 16]. 

4.2 The rescaled equations 

In Sec. 3.3 we showed that, to linear order the peculiar motion of the matter is also described 
by a rescaled velocity field which is both acceleration-free and irrotational. The introduction of 
the rescaled peculiar velocity Va = a'^^'^Va, leads immediately to the rescaling of the irreducible 
kinematical variables of the peculiar motion. In particular, given that t)bVa = o}^^i)bVa, we have 

e = a^/^B , dab = a^^^^ab and LOab = a'^'^Qah , (36) 

where 

Q = taV\ c7a6 = D(5K) and cUa6 = D[6K] (37) 

are respectively the rescaled peculiar contraction scalar, the recsaled peculiar shear and the 
rescaled peculiar vorticity. 

On using expressions (36) one can transform Eqs. (32)- (35) into the following system 

V(a> = -a-^/2^(,) - iey^,) , (38) 
Oah = -|a^/^ea„5-a^/^f^c(af^6>'-a^/^<^c(a'^6)'' + D(„H), (40) 

^ab = -la^/^eUab + 2a^^^dc[a^bf-%Vb], (41) 

which governs the motion of the matter, as described by the rescaled peculiar velocity field Va, 
to second order. 

5 The Zeldovich approximation 
5.1 Propagation equations and constraints 

The Zeldovich approximation addresses the mildly nonlinear collapse of protogalactic clouds, 
on scales well within the Hubble radius, as they decouple from the background expansion and 
turn around. The approximation builds upon the exact linear result of acceleration-free and 
irrotational motion of the dust component (see Sec. 3.3), which is extrapolated into the nonlinear 
regime. This considerably simplifies the equations and allows one to obtain analytical solutions. 
Here, we will describe the mildly nonlinear regime through the set of the rescaled second order 
equations (38)- (41). When the Zeldovich ansatz Va = = Uah is applied, the motion of the 
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collapsing pressure-free matter is governed by Eqs. (39) and (40) only, which for constant velocity 
and in the absence of vorticity reduce even further to 




(42) 
(43) 



At the same time, Va = = coab means that Eq. (38) transforms into the constraint 



(44) 



and that Eq. (41) becomes an identity. The latter ensures that setting the vorticity to zero at 
second order is a self-consistent assumption. This is not the case, however, with constraint (44) 
whose consistency in time is not straightforward. Nevertheless, constraint (44) does not pose any 

real problem for the consistency of the Zeldovich ansatz at second order. Indeed, by propagating 
(44) in time one obtains the following evolution equation for the 4-acceleration vector 



Since the zero-vorticity constraint is automatically satisfied at second order, the above propa- 
gation equation is effectively the consistency condition of the Zeldovich approximation to this 
perturbative order. 

We now turn our attention to the matter term at the end of Eq. (42) , which conveys the effects 
of spacetime curvature (see Sec. 2.3). Given that Va is constant and that /x oc a~^, the impact 
of the background matter upon decays away. In other words, as the collapse is progressively 
dominated by the kinematics and the role of gravity becomes negligible. The situation is closely 
analogous to that observed in studies of silent universe models (e.g. see [16]), or in Bianchi I 
cosmologies, where the latter become shear dominated at early times (e.g. see [19]). 

One should note that the effect of the approximation has been to decouple the evolution along 
neighboring world lines from each other, similar to what happens in silent universe models, 
where (in an inhomogeneous situation) one can integrate the evolution equations separately 
along neighboring worldlines. In other words, the system of partial differential equations has 
been reduced to an effective system of ordinary differential equations. 

5.2 The shear eigenframe 

To proceed further we introduce the new "time" variable r, constructed so that f = —a 
where the minus sign compensates for the fact that we are dealing with a collapsing region 
(i.e. < 0). This choice guarantees that f > always and that r ^ oo as we approach 
the singularity and — — oo. Then, Eqs. (42), without matter, and (43) transform into the 
following set 



where from now on a prime will indicate differentiation with respect to r. Assuming the shear 
eigenframe we have aab = (o'u, 0^22 , o'ss), where (J33 = — ((Tii+(722) due to the trace-free nature of 
the shear. Then, if ai = an, 02 = 012 and aa = 0^33 determine the three shear eigen-directions, 



\QAa + \o^l'^^lVa . 



(45) 



e' = |e + 29-^^2 




(46) 
(47) 
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the second order system (46), (47) reads 



e' = l@ + 2e~\af + a1 + aia2) , (48) 
<fi = |cTi + ie-ifj?-|e-i(ai + CT2)a2, (49) 

= 1^2 + |0"^^i - |e~H^i + ^2)^1 , (50) 

where the behavior of CJ3 is determined from that of ai, a2 entirely. The above second-order 
equations govern the small-scale evolution of the pressure-free matter, as the latter decouples 
from the background expansion and starts to turn around, and provide a fully covariant formu- 
lation of the Zeldovich approximation. 

The adoption of the shear eigenframe means that the three shear eigen-directions acquire 
special importance. Consequently, it is convenient to introduce the quantities 

ei = ai + lQ, (51) 

with i = 1, 2, 3, which provide the contraction rates in each shear eigen-direction. They will also 
be used to define a scaling length along each direction (see Sec. 5.4). 



5.3 The attractors 

The question is whether or not the relativistic analysis also predicts that one-dimensional pan- 
cakes are the final configurations of any generic collapsing ovcrdcnsity. Given the qualitative 
nature of the question one can employ a dynamical system approach to provide the answer. We 
begin by defining the following dimensionless dynamical variables (see also [27, 28]) 

^^^3(2L±M and E^ = ^^l(5glM, (52) 

which in turn imply the auxiliary relations 

ai = ie(E+-F V^E_), ^2 = 5e(E+ - V3E_) and a3 = -|eE+. (53) 

Clearly, E+ and E_ measure the kinematical anisotropy of the collapse. If both are zero the 
collapse is spherically symmetric. When only E_ = we have the degenerate case case of two 
equal shear eigenvalues. On introducing T,± Eq. (48) transforms into 

e' = ie + §e (e^ + e^ ) , (54) 

while Eqs. (49) and (50) read [27] 

EV = i [1 - E+ - 2(E^ + E^ )] E+ + iE?. , (55) 
E'_ = i [1 + 2E+ - 2(e2 + E?_)] E_ . (56) 

Accordingly, the evolution of E+ and E_ has decoupled from that of and the shape evolution 
of the overdensity is now monitored by the sub-system (55), (56). Technically speaking, the 
problem has now been reduced to the study of the planar dynamical system (55), (56) depicted 
in Fig. 1. Physically, the dimensional reduction means that the shape of the collapsing dust 
cloud does not depend on the time scale of the collapse. Note that our equations are invariant 
under the cyclic symmetry 
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Figure 1: Phase plane S+ = X, S_ = Y. The Hnes forming the central triangle correspond 
to ai = —0/3 (with i = 1,2,3). The three pancakes arc located at the intersections of these 
lines. The lines bisecting the vertices correspond to two equal shear eigenvalues, (e.g. E_ = 
corresponds to ai = [27]. 



(Ti ^ (72 ^ 0-3 ^ cri , (57) 

which is equivalent to the chain 

— ^ — — — ^ — "I" — ^ ^+ ' (58) 
S_ ^ +^S+-iE_ ^ -^S+-iE_ ^ S_. (59) 

The two sets of lines in Fig. 1, those that form the central triangle and those that bisect 
the vertex angles, are each invariant submanifolds (mapped into themselves by the flow). The 
bisecting lines represent the rotationally symmetric solutions with trajectories 

LI: S_ = ^ ai = a2, (60) 
L2: S_ = V3S+ <^ ^2 = ^3, (61) 
L3 : S_ = - V3S+ <^ ai = a3. (62) 

They are mapped into themselves by the symmetries noted above, so they are equivalent to each 
other. Concentrating on LI, we notice that it is involutive because S_ = immediately implies 
S'_ = (see Eq. (56)). The equation along this trajectory is 

= i [1 - S+ - 2E^] S+ , S_ = (63) 

representing a family of rotationally symmetric solutions. It has three fixed (stationary) points, 
which are each self-similar solutions of the system (55), (56). In particular we have 

PI (S_ = 0, E_|_ = —1): In this case ai = a2 = —0/3 and = 26/3, which implies that 
©1 = ©2 = and ©3 = 0. This point corresponds to a rotationally symmetric pancake collapse. 
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It occurs at a triangle vertex, where two of the triangle generating lines and a bisecting line all 
intersect, and is an attr actor in this orbit. 

P2 (E_ = 0, E+ = i): Here 0-1 = ^2 = G/6 and CT3 = -6/3 meaning that Gi = 62 = 6/2 
and 03 = 0. This point corresponds to a spindle or filamentary rotationally symmetric collapse. 
It occurs where a bisecting line meets a triangle generating line, and is an attractor in the orbit. 

P3 (E_ = 0,E+ = 0): At this point ai = ^2 = CT3 = and Gi = 62 = 63 = 0/3, which 
corresponds to an isotropic collapse. It occurs at the centre where the three bisecting lines meet, 
and is a source in the orbit. Note that the same set of fixed points also occur on each of the 
other two bisecting lines L2 and L3. 

The triangle generating lines represent the (generally anisotropic) 2-dimensional collapse 
solutions with trajectories 

LA: 0i = cTi + i0 = O 4» E_ = -^(E+ + 1), (64) 
L5 : 02 = (72 + i© = ^ ^- = +71 (^+ + 1) ' (65) 
L6 : 03 = aa + 5© = O ^+ = \ ■ (66) 

They are mapped into themselves by the symmetries noted above, so they are equivalent to 
each other. Let us concentrate on L6, which is involutive since E+ = 1/2 leads to E^ = (see 
Eq. (55)). The equation along this trajectory is 

S+ = I> 5]'_ = |[f-E2_]E_, (67) 

representing a family of exact 2-dimensional collapse solutions. The latter are rotationally 
symmetric only where they intersect the bisecting lines. This orbit has three fixed points. The 
point (E_ = 0, E+ = 1/2), where it intersects a bisecting line, namely the LRS spindle P2 noted 
above. Also, the two points (E_ = '\/3/2, E_|_ = 1/2) and (E_ = — \/3/2, E_|_ = 1/2), where L6 
meets another generating line and a bisecting line simultaneously. Each of the latter two points 
is equivalent to the pancake collapse node PI noted above. 

In summary, the vertices of the triangle are stationary points of the system (55), (56) and 
are attractors. Each is a one-dimensional pancake solution, with a stationary state along two 
directions and collapse along the third shear eigen-direction. Generic solutions are asymptotic to 
one of these three one-dimensional pancakes (one for each eigen-direction) and so are asymptotic 
to one-dimensional pancakes. The bisecting lines intersect at the point that corresponds to a 
shear-free spherically symmetric collapse, also a stationary point. Where the bisecting lines 
intersect the triangle, we have stationary points that represent exact filamentary solutions. 
The pancakes are stable nodes, the filaments are saddle points, and the spherically symmetric 
collapse is an unstable node. This proves that, once collapse sets in, the pancakes are the natural 
attractors for a generic overdensity. 

5.4 The stationary points 

In cosmology we use the average expansion scalar to define the cosmological scale factor. Simi- 
larly, one can utilize the collapse rate (0) of an inhomogeneity to define an average local scaling 
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length {£) as well as directional scaling lengths (4, i = 1,2,3) along the three shear eigen- 
directions (e.g. see [16]). In particular we have 

j = l@ and ^ = ai + le = @i. (68) 

Thus, £ is the geometric mean of the three directional scale factors ii, while @ is the average 
of their expansion rates. Following the standard terminology of [29], pancakes occur when two 
of the iiS approach finite vales while the other goes to zero, spindle-like singularities correspond 
to two of the iiS going to zero and the third to a constant, and finally a point-like singularity 
occurs when all three of the iiS tend to zero. 



As we have already pointed out, stationary points correspond to finite values of the dimen- 
sionless dynamical variables T,± introduced in Sec. 5.3. Concentrating on the three stationary 
points along the LI trajectory, we can calculate the associated collapse timescales. In particular, 
following closely the analysis given in [16] we have: 



PI (S_ = 0, = —1): This case corresponds to pancake collapse along the third shear 
eigen-direction, Eq. (54) implies that O' = G. Here, the evolution of the contraction scalar G is 
given by 

e = . (69) 

where t is the proper time along the matter frame and the zero suffix indicates the beginning 
of the collapse. Note that we have used the transformation law (36a) and assumed a dust- 
dominated universe (i.e. a oc t"^^^). Then, by integrating Eqs. (68) we find that the average scale 
factor and the scaling length along the direction of the collapse evolve as 



1 + (^4/3 _ ^4/3) - - ^ ^ (^4/3 _ 



1/3 



(70) 



respectively. Note that ii, I2 = constant and ai = S@o/4t^/^ < 0, given that ©0 < for 
collapse. Thus, the pancake singularity occurs at 



.1/4 



{to - |©o 



3/4 



(71) 



P2 (I]_ = 0, = ^): For filamentary collapse along the first two shear eigen-directions 
Eq. (54) gives 9' = G/2. This leads to the following evolution law 



e = 



8Gnt^/^ 



84/' + 3Go(t4/3-t^^ 



4/3 X 



(72) 



for the average contraction scalar. In addition, we find that the average and the two directional 
scaling lengths evolve as 



h = {h)o l + a2(i'/'-iJ/') 



and io 



2/3 



^^2 



1 + 02(^^/3 -if) 



(73) 
(74) 
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respectively, while £3 = constant (0:2 = 30o/8i^/^). According to the above, the spindle-like 
singularity takes place at 

U=t'J'{to-l%'y^\ (75) 



P3 (S_ 
implies 



0,Eh 



G 



0): Finally, for isotropic collapse 8' 



4/3 X 



and £ = £1 



= 0/3 (see Eq. (54)) which in turn 



for the contraction scalar and the scaling length respectively [a^ 
singularity occurs at 

, _,l/4/. .a_A3/4 



(*o-4eo-^)^ 



(76) 



9o/4i^/3). Here, the point 



(77) 



Clearly, the contracting behavior at each stationary point is monotonic and fixed by the 
initial conditions. Given that 60 < 0, there is always a singularity in the future. According to 
Eqs. (71), (75) and (77), for the same set of initial conditions, the pancake singularity occurs 
first and the point singularity last. This result, which is due to the increase in the contraction 
rate of fluid congruence in the presence of shear, is in agreement with the qualitative "collapse 
theorem" given in [15]. 



6 Conclusions 

A key issue in contemporary theoretical cosmology is understanding the physics of gravitational 
collapse and the mechanisms that have given rise to the observed large-scale structure of the 
universe. Linear perturbations, about a Friedmann-Robertson- Walker model, are fairly straight- 
forward to follow. When certain simplifying symmetries are imposed, the nonlinear collapse can 
also be treated analytically. The spherical top-hat model, in Newtonian theory, and the Tolman- 
Bondi solution, in General Relativity, are such examples. When it comes to the more realistic 
non-spherical collapse, however, the Zeldovich approximation has been the most influential and 
celebrated paradigm. It has also inspired a number of variations, all of which try to address the 
complexities of structure evolution beyond "caustic" formation. Despite the number of these 
different treatments, however, a fully covariant approach to the Zeldovich approximation has 
been missing. This is the issue that the present paper has tried to address. 

Wc have pursued a relativistic treatment of nonlinear peculiar velocities by adopting a two- 
frame approach which enables us to provide a truly Lagrangian formalism. Assuming matter 
with non-relativistic peculiar motion relative to a quasi-Newtonian frame, we have derived the 
fully nonlinear equations and applied them to the case of pressure- free dust. The inherent 
transparency of the covariant formalism means that all the terms in our equations have a clear 
physical and geometrical interpretation. For a given background, it is also straightforward to 
identify the perturbative order of all the effects that influence the peculiar motion of the matter. 
Assuming an Einstein-de Sitter background and allowing for up to second order perturbative 
terms, we have addressed the mildly nonlinear collapse and in the process provided a fully 
covariant formulation of the Zeldovich approximation. On introducing the Zeldovich ansatz 
of acceleration-free and irrotational peculiar motion, our equations have reduced to a planar 
dynamical system, which has one-dimensional pancakes as its natural attractors. Thus, just like 
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the Newtonian treatment, the relativistic approach also predicts that, as long as we consider the 
mildly nonlinear stage, pancakes are the final end-states of any generic collapsing overdensity. 
In addition, by looking at the collapse timescales associated with the stationary points, we have 
found that pancake collapse takes place first and point-like singularities occur last. 

Our results for the fate of generic collapsing overdensities are in disagreement with stud- 
ies of silent-universe dynamics, which argue for spindle-like rather than pancake singularities 
(see [15, 16]). The reason for this difference lies in the role of the tidal field within the adopted 
cosmological model. Silent universes allow for non-zero electric Weyl tensor, but set its mag- 
netic counterpart to zero. In the presence of this "truncated" tidal field the fluid element evolves 
towards a Kasner-type singularity, where pancakes are a set of measure zero. Our two-frame 
analysis of a perturbed Einstein-de Sitter model, on the other hand, shows that the influence of 
the long-range gravitational field is negligible at second order. Technically speaking, the covari- 
ant equations governing the mildly nonlinear collapse of a dust cloud have zero Weyl-curvature 
input, effectively reducing to the Newtonian ones. As a result, pancakes are reinstated as the 
natural attractors of a generic collapsing overdensity, this time within the relativistic framework. 
This is in agreement with numerical simulations which clearly favor pancake formation over all 
other types of structures (see [14]). Once again we point out that our result holds at second 
perturbative order and refers to what is known as the mildly nonlinear regime. Although the 
second order equations can probe deep into the epoch of structure formation, they are expected 
to break down close to the singularity. Clearly, if we were to address the full picture, we need 
to incorporate the tide effects. For an Einstein-de Sitter background, however, the electric and 
of the magnetic Weyl terms are both of third perturbative order (see Eq. (22)). In other words, 
the fully nonlinear collapse evolves under the influence of the full tidal field rather than the 
purely-electric Weyl field of the silent models. In addition, at higher than the second perturba- 
tive order, the magnetic Weyl tensor is also a source of peculiar vorticity (see Eqs. (18), (23)) 
and one can no longer consistently argue for irrotational peculiar motion. Of course, as the 
collapse proceeds beyond the mildly nonlinear stage, pressure gradients also become important, 
which means that an acceleration-free motion is no longer sustainable, and the whole idea of the 
Zeldovich approximation is expected to break down. 

Note that, although it is usually bypassed, vorticity is an important issue, given the ob- 
served rotation of galaxies and galaxy clusters. By construction, the Zeldovich approximation 
cannot address this question. Newtonian modifications of the standard approach to incorporate 
vorticity have already been suggested in the past (see [8]). The covariant formalism developed 
here can provide the basic framework for a relativistic, mathematically rigorous and physically 
transparent treatment of rotating dust during the nonlinear collapse. 
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